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SHEAR-FLEXIBLE  FINITE-ELEMENT  MODELS  OF  LAMINATED 
COMPOSITE  PLATES  AND  SHELLS 

Ahmed  K.  Noor*  and  Michael  D.  Mathers* 

Langley  Research  Center 

SUMMARY 

Several  finite-element  models  are  applied  to  the  linear  static,  stability,  and  vibration  analy¬ 
sis  of  laminated  composite  plates  and  shells.  The  study  is  based  on  linear  shallow-shell  theory, 
with  the  effects  of  shear  deformation,  anisotropic  material  behavior,  and  bending-extensional 
coupling  included.  Both  stiffness  (displacement)  and  mixed  finite-element  models  are  considered. 
Discussion  is  focused  on  the  effects  of  shear  deformation  and  anisotropic  material  behavior  on 
the  accuracy  and  convergence  of  different  finite-element  models.  Numerical  studies  are  pre¬ 
sented  which  show  the  effects  of  (a)  increasing  the  order  of  the  approximating  polynomials, 

(b)  adding  internal  degrees  of  freedom,  and  (c)  using  derivatives  of  generalized  displacements 
as  nodal  parameters. 


INTRODUCTION 

Although  the  finite-element  analysis  of  isotropic  plates  and  shells  has  received  considerable 
attention  in  the  literature,  investigations  of  laminated  composite  plates  and  shells  are  rather 
limited  in  extent.  The  reliable  prediction  of  the  response  characteristics  of  high-modulus 
fibrous  composite  plates  and  shells  often  requires  inclusion  of  the  transverse  shear  effects  in 
their  mathematical  models.  This  fact  has  been  amply  documented  for  linear  static,  stability, 
and  dynamic  problems.  (See,  for  example,  refs.  1  to  5.) 

At  present  there  are  three  approaches  for  developing  plate  and  shell  finite-element  models 
which  account  for  shear  deformation.  The  first  approach  is  based  on  the  use  of  three- 
dimensional  isoparametric  solid  elements  which  automatically  include  the  shear-distortion  mecha¬ 
nism  (refs.  6  and  7).  The  second  approach  employs  two-dimensional  elements  used  with  inde¬ 
pendent  shape  (or  interpolation)  functions  for  displacements  and  rotations  (refs.  8  and  9).  The 
third  approach  is  based  on  the  addition  of  effects  of  shear  deformation  to  two-dimensional 
classical  plate  or  shell  elements  through  the  use  of  equilibrium  equations  (refs.  10  and  11). 
Although  it  is  desirable  to  have  an  element  which  gives  accurate  results  regardless  of  how 
important  the  shear  deformation  is,  most  of  the  existing  elements  do  not  satisfy  this 
requirement. 
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In  the  context  of  the  stiffness  method,  the  first  approach  has  the  major  disadvantage 
that  it  leads  to  a  stiffness  matrix  which  is  (1)  very  large  for  laminated  composites  consisting 
of  many  layers  and  (2)  highly  ill  conditioned  for  thin  plates  or  shells.  If  low-order  interpola¬ 
tion  polynomials  are  used,  the  second  approach  leads  to  overly  stiff  elements  for  very  thin 
plates  and  shells.  Although  the  aforementioned  drawbacks  have  been  recognized  and  some 
improvements  have  been  suggested,  the  difficulties  have  not  been  overcome.  (See,  e.g., 
refs.  12  to  17.)  The  range  of  validity  of  the  third  approach  has  not  been  explored.  Since 
the  second  approach  provides  flexibility  and  simplicity  in  fulfilling  the  interelement  compati¬ 
bility  conditions  and  does  not  result  in  as  large  a  stiffness  matrix  as  in  the  first  approach,  it 
was  adopted  in  the  present  study. 

The  first  objective  of  this  paper  is  to  assess  the  relative  merits  of  a  number  of  displace¬ 
ment  and  mixed  shear-flexible  finite  elements  when  applied  to  the  linear  static,  stability,  and 
vibration  problems  of  laminated  plates  and  shells.  Emphasis  is  focused  on  the  effects  of  shear 
deformation  and  anisotropic  material  behavior  on  the  accuracy  and  convergence  of  the  different 
models.  The  second  objective  is  to  study  the  effects  of  increasing  the  order  of  approximating 
polynomials,  adding  internal  degrees  of  freedom,  and  using  derivatives  of  generalized  displace¬ 
ments  as  nodal  parameters  on  the  accuracy  and  rate  of  convergence  of  the  different  models. 

To  the  authors’  knowledge  no  publication  exists  in  which  the  aforementioned  effects  are 
studied  in  any  detail. 

The  analytical  formulation  is  based  on  a  form  of  the  shallow-shell  theory  modified  to 
include  the  effects  of  shear  deformation  and  rotary  inertia.  Indicial  notation  is  used  through¬ 
out  this  paper  since  it  is  particularly  useful  in  identifying  the  symmetries  and,  consequently, 
simplifies  the  element  development.  Both  triangular  and  quadrilateral  elements  are  considered. 
The  elements  are  conforming  and  satisfy  continuity  requirements  of  the  type  (continuity 
of  the  fundamental  unknowns). 


SYMBOLS  AND  NOTATION 


shell  compliance  coefficients,  inverse  of  shell  stiffnesses 

^a^yp^^aj^jp  J 

a  side  length  of  plate  or  shallow  shell 

^aPyp  extensional  stiffnesses  of  shell 

transverse  shear  stiffnesses  of  shell 


(k)  (k) 

^aPyp’^aSPS 


stiffness  coefficients  of  kth  layer  of  shell 
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consistent  nodal  load  coefficients 
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W,W  work  done  by  external  forces 

X3  orthogonal  curvilinear  coordinate  system  (see  fig.  1) 

xj^  nodal  values  of  x^ 

^  dimensionless  eigenvalues  of  stiffness  matrix 

fj.  relative  size  of  rth  element  in  variable  grid  (eq.  (37)) 

?7r,i?r+i  dummy  coordinates  of  ends  of  rth  element 

d  fiber  orientation,  angle  between  fiber  direction  and  xj-axis 

K  constant  defined  in  appendix  D 

X  in-plane  loading  parameter 

X  nondimensional  frequency  ^co\/pa^y^Ej  for  plates;  for  shallow 

spherical  segments;  cu  \  ff)^ for  circular  cylinders^ 

V  Poisson’s  ratio  for  isotropic  materials 

^LT  Poisson’s  ratio  measuring  strain  in  T-direction  (transverse)  due  to  uniaxial  normal 

stress  in  L-direction  (direction  of  fibers) 

natural  coordinates  of  node  i 


n,nj^ 


natural  (dimensionless)  coordinate  system  in  element  domain 

functionals  defined  in  equations  (1)  and  (2) 

density  of  plate  or  shell  material 

density  of  kth  layer  of  laminated  shell 

uniform  extensional  stress  in  cylindrical  shell 


rotation  components 
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nodal  displacement  parameters 


shell  domain 


CO  circular  frequency  of  vibration  of  shell 


a: 


Range  of  indices: 

Lowercase  Latin  indices  1  to  m 

Uppercase  Latin  indices: 

I,J  1  to  5 
I,J  1  to  8 

Greek  indices  1 ,2 

Finite-element-model  notation: 

SQN  stiffness  formulation,  quadrilateral  element,  N  shape  functions  per  fundamental 

unknown 

STN  stiffness  formulation,  triangular  element,  N  shape  functions  per  fundamental 

unknown 

MQN  mixed  formulation,  quadrilateral  element,  N  shape  functions  per  fundamental 

unknown 

MTN  mixed  formulation,  triangular  element,  N  shape  functions  per  fundamental  unknown 

SQH  stiffness  formulation,  quadrilateral  element,  Hermitian  interpolation  functions 

MATHEMATICAL  FORMULATION 

The  analytical  formulation  is  based  on  a  form  of  the  shallow-shell  theory,  with  the 
effects  of  shear  deformation,  anisotropic  material  behavior,  rotary  inertia,  and  bending- 
extensional  coupling  included.  (See  appendix  A  and  ref.  18.)  For  stability  problems,  the 
prebuckling  stresses  are  assumed  to  be  given  by  the  momentless  (membrane)  theory.  Two 
finite-element  formulations  are  considered.  In  the  first  formulation  (displacement  model)  the 
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fundamental  unknowns  consist  of  the  displacement  and  rotation  components  of  the  shell 
reference  (middle)  surface,  and  the  stiffness  matrix  is  obtained  by  using  Hamilton’s  principle 
(which  for  static  problems  reduces  to  the  principle  of  minimum  potential  energy).  The  funda¬ 
mental  unknowns  in  the  second  formulation  (mixed  model)  consist  of  the  13  shell  quantities: 
generalized  displacements  u^,,  w,  and  and  stress  resultants  and  (See 

fig.  1  for  sign  convention.)  The  generalized  stiffness  matrix  is  obtained  by  using  a  modified 
form  of  the  Hellinger-Reissner  mixed  variational  principle. 

The  functionals  used  in  the  development  of  displacement  and  mixed  models  are  given  by 
the  following  equations: 


Displacement  models 

-  U  +  -  W  -  T  (1) 

Mixed  models 

=  V  +  U°  -  -  W  -  W  -  T  (2) 

where 


U  -  y  J" {CQ;j3.yp  (dQ;U|3  9-yUp  +  2k^p  ^  ^^a^yp(pa^^  ^y^P 

+  k^^  d^<l>pV/)  +  D^^p  a^0p  +  Cc,303(aaW  dpw  +  2<t>a  d^w  +  (3) 

U°  =  ^  9q,w  dpw  dfl  (4) 

V  [Na/)(9aU^  +  w)  +  d^<pp  +  Q^(<t>o,  +  aaw)]df2  (5) 

U*“  -  ^  N.yp  +  2Bq^Pjp  NQ,|g  ^7P  ^oc^yp  ^7P  ■^0(3j33  ^oc  Q/3}4f2  (6) 

^  ^  2’^'^  J  [^o(%^a;  2mju^0Q,  +  m20Q,<^^4S2  (7) 

W  "/{Pcx^cx  +  Pw)df2  +  +  Qpw)npdc  (8) 

W  +  %)  +  ^ap(-'^a  +  ^a)  +  +  wjjn^dc  (9) 


In  equations  (3)  to  (9),  C^pyp,  ^a^yp’  F^^.^p  are  extensional  stiffnesses,  bend¬ 
ing  stiffnesses,  and  stiffness  interaction  coefficients  of  the  shell;  traiisverse  shear 
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stiffnesses  of  the  shell;  ^aPyp^  ^aj5jp^  ^a3(53  compliance  coefficients 

(see  appendix  A);  are  the  curvature  components  and  twist  of  the  shell  surface;  XN^l3 

are  the  initial  stress  resultants  (prestress  field)  which  are  proportional  to  the  in-plane  load  fac¬ 
tor  X;  p^  and  p  are  the  external  load  components  in  the  orthogonal  coordinate  direc¬ 
tions  x^  and  X3,  respectively;  niQ,  mj,  and  are  density  parameters  of  the  shell 

defined  in  appendix  B;  co  is  the  circular  frequency  of  vibration  of  the  shell;  SI  is  the 
shell  domain;  c^  and  c^^  are  portions  of  the  boundary  over  which  tractions  and  displace¬ 
ments  are  prescribed;  n^  is  the  unit  normal  to  the  boundary;  the  quantities  with  a  tilde 

denote  prescribed  boundary  stress  resultants  and  displacements;  and  3^  =  — . 

dx^ 

FINITE-ELEMENT  DISCRETIZATION 

The  shell  region  is  decomposed  into  finite  elements  connected  at  appropriate 

nodes,  where  the  superscript  e  refers  to  the  element.  A  typical  element  is  isolated  from  the 
model  and  the  fundamental  unknowns  are  approximated  by  expressions  of  the  form: 


Displacement  models 
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(10) 

CO 

II 

(11) 

(12) 

Mixed  models 

In  addition  to  the  approximations  of  the  generalized  displacements  (eqs. 
stress  resultants  are  approximated  by 

(10)  to  (12)),  the 

Nai3  = 

(13) 

(14) 

(15) 

where  superscripts  identify  the  location  and  subscripts  designate  the  ordering  of  nodal  unknowns; 

are  the  shape  (or  interpolation)  functions;  \I/j  (i  =  1  to  m,  J  =  1  to  5)  are  nodal  dis¬ 
placement  parameters  (including,  possibly,  nodeless  variables);  Hj  (i  =  1  to  m,  J  =  1  to  8) 
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are  nodal  stress-resultant  parameters;  m  equals  the  number  of  shape  functions  in  the  approxi¬ 
mation;  Greek  indices  take  the  values  1,2;  and  a  repeated  lowercase  Latin  index  denotes  sum¬ 
mation  over  the  range  1  to  m. 

ELEMENT-BEHAVIOR  REPRESENTATION 

A  number  of  displacement  and  mixed  finite  elements  having  both  triangular  and  quadri¬ 
lateral  shapes  were  developed  in  the  present  study.  All  the  elements  satisfy  the  continuity 
conditions  required  by  the  variational  principles  on  which  they  are  based.  Within  each  family 
of  elements,  different  shape  (or  interpolation)  functions  are  used  for  approximating  the  funda¬ 
mental  unknowns.  The  characteristics  and  designations  of  these  elements  are  summarized  in 
table  1  and  are  referred  to  frequently  in  the  subsequent  sections. 

All  the  triangular  elements  developed  are  based  on  complete  polynomial  approximations 
of  the  fundamental  unknowns,  thus  ensuring  that  the  functional  variation  is  independent  of 
coordinate  transformations.  Most  of  the  quadrilateral  elements  considered  in  the  present  study 
are  of  the  serendipity  type  (refs.  19  and  20),  that  is,  with  their  nodes  located  along  the  ele¬ 
ment  boundaries.  The  polynomial  approximations  used  in  these  elements  include  terms  which 
are  of  higher  order  than  the  complete  expansion,  and  therefore,  the  functional  variation  is 
dependent  on  coordinate  transformation. 

In  each  element,  the  same  set  of  shape  functions  is  used  for  approximating  all  the  fun¬ 
damental  unknowns  and  the  nodal  parameters  are  selected  to  be  the  values  of  the  fundamen¬ 
tal  unknowns  at  the  different  nodes.  However,  in  one  of  the  elements  (SQ8-4  element), 
polynomials  of  different  degree  were  used  for  approximating  different  sets  of  fundamental 
unknowns  (lower  degree  polynomials  were  used  for  approximating  the  rotations);  in  the 
SQH  element,  products  of  first-order  Hermitian  polynomials  were  chosen  as  shape  functions 
and  the  nodal  parameters  consisted  of  the  generalized  displacements,  their  first  derivatives,  and 
mixed  second  derivative  with  respect  to  the  dimensionless  local  coordinates  and  ^2- 

(See  appendix  C.)  Continuity  of  these  derivatives  is  enforced  along  the  interelement  bound¬ 
aries.  Since  this  is  not  required  by  the  variational  principle,  the  element  is  overconforming. 

For  the  two  quadrilateral  stiffness  elements  with  four  and  eight  nodes,  internal  degrees 
of  freedom  are  added  through  the  addition  of  displacement  modes  which  vanish  along  the 
edges  of  the  element.  Those  modes  are  usually  called  bubble  functions  (ref.  21).  The  shape 
functions  associated  with  the  internal  degrees  of  freedom  are  products  of  the  equations  of  the 
element  boundaries  times  another  polynomial,  with  the  product  representing  bubble  or  internal 
displacement  modes  (elements  SQ5,  SQ7,  SQ9,  and  SQll).  The  case  of  one  internal  mode 
(SQ5  and  SQ9)  corresponds  to  zero  degree  of  the  latter  polynomial.  (See  table  1  and 
appendix  C.) 
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In  all  the  elements  developed,  the  rigid  body  modes  that  cause  no  straining  have  not 
been  included  explicitly  in  the  displacement  fields;  rather,  implicit  representation  of  these  modes 
was  made.  A  quantitative  estimate  of  the  accuracy  of  rigid-body-mode  representation  was  made 
by  evaluating  the  six  lowest  eigenvalues  of  the  element  stiffness  matrix.  This  is  discussed  fur¬ 
ther  in  connection  with  the  numerical  studies. 

For  modeling  shells  with  curved  boundaries,  isoparametric  elements  were  used  in  which 
the  element  boundary  curves  are  approximated  by  the  same  shape  functions  used  in  approxi¬ 
mating  the  behavior  functions,  that  is, 

(16) 

where  x^  are  the  nodal  values  of  x^.  Numerical  results  obtained  with  the  use  of  isopara¬ 
metric  SQ12  elements  are  presented  in  the  next  main  section. 

FINITE-ELEMENT  EQUATIONS 


The  governing  equations  for  each  element  are  obtained  by  first  replacing  the  fundamental 
unknowns  by  their  expressions  in  terms  of  the  shape  functions  (eqs.  (10)  to  (15))  in  the 
appropriate  functional  (action  integral  for  displacement  models  and  Hellinger-Reissner  functional 
for  mixed  models)  and  then  applying  the  stationary  conditions  of  that  functional.  This  leads 
to  a  set  of  equations  for  each  element  of  the  following  form: 


Displacement  models 

K«J  +  XKf]  =  p;  + 

V-  - ' 

Mixed  models 


(17) 


where  and  are  stiffness  and  geometric,  or  initial  stress,  stiffness  coefficients; 

ii  IJ  IJ  ii  -ii 

M  are  consistent  mass  coefficients;  S—  and  S-  are  “generalized”  stiffness  coefficients; 

IJ  .  IJ  IJ 

and  Pj  are  consistent  load  coefficients.  Tlie  formulas  for  tlie  aforementioned  stiffness,  mass, 
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and  load  coefficients  are  given  in  appendix  D.  For  stress-analysis  problems,  X  =  co  -  0;  for 
bifurcation-buckling  problems,  co  =  Pj  =  0;  and  for  free-vibration  problems,  X  =  Pj  =  0. 

In  equations  (17)  and  (18)  the  range  of  the  lowercase  Latin  superscripts  is  1  to  m;  the 
range  of  the  uppercase  Latin  subscripts  (I,J)  and  (i,J)  is  1  to  5  and  1  to  8,  respectively. 

The  K,  M,  and  S  terms  are  completely  symmetric  under  the  interchange  of  one  pair  of 
indices  for  another,  each  pair  of  indices  consisting  of  a  superscript  and  a  subscript  just  beneath 
it. 

To  write  equations  (17)  and  (18)  in  matrix  form,  the  first  superscript-subscript  pair  of 
each  of  the  K,  S,  and  M  terms  defines  the  row  number  and  the  second  pair  defines 
the  column  number.  For  example,  in  equations  (17)  the  term  kJj  is  located  in  the 
[5(i-l)  +  ijth  row  and  the  fSCj-l)  +  J]th  column  of  the  element  stiffness  matrix. 

In  the  stress-analysis  problems,  the  internal  degrees  of  freedom  (nodal  parameters  associ¬ 
ated  with  bubble  modes)  can  be  eliminated  without  any  loss  of  accuracy  by  using  the  static 
condensation  procedure  (ref.  22).  In  stability  and  vibration  problems,  this  is  not  done  since 
it  results  in  approximate  elemental  matrices. 

The  integrals  in  the  expressions  for  the  stiffness,  mass,  and  load  coefficients  (appendix  D) 
are  evaluated  by  means  of  the  numerical  quadrature  formulas  presented  in  references  20  and  23. 
In  each  case,  the  quadrature  formula  selected  had  the  least  number  of  points  required  to  ensure 
exact  evaluation  of  the  integrals  (depending  on  the  degree  of  the  interpolation  polynomials). 
Exceptions  to  this  are  the  cases  of  general  quadrilateral  or  isoparametric  elements  based  on  the 
displacement  models  in  which  the  stiffness  and  geometric  stiffness  coefficients  contain  fractional 
rational  functions  that  are  approximated  by  polynomials  in  the  numerical  quadrature  process. 
Each  entry  in  the  elemental  matrices  S  and  S  of  the  mixed  models  (eqs.  (18))  contains 
just  a  single  term.  (See  appendix  D.)  In  contrast,  the  entries  of  the  matrix  K  of  the 
displacement  models  (eqs.  (17))  are  linear  combinations  of  at  least  four  terms,  as  implied  by 
the  repeated  (dummy)  subscripts  of  the  coefficients  K  in  appendix  D.  In  view  of  this,  the 
formation  of  the  elemental  matrices  for  the  mixed  models  is  simpler  and  was  found  to  be  less 
time  consuming  than  for  the  displacement  models. 

BOUNDARY  CONDITIONS 

In  the  displacement  models,  only  kinematic  (geometric)  boundary  conditions  need  to  be 
satisfied.  Eorce  (stress)  boundary  conditions  can  also  be  satisfied  if  displacement  derivatives 
are  chosen  as  nodal  parameters  (e.g.,  SQH  element).  The  effect  of  introducing  the  stress 
boundary  conditions  on  the  accuracy  of  solutions  is  discussed  in  the  examples  in  the  section 
“Numerical  Studies.” 
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In  the  mixed  models,  both  kinematic  and  force  (stress)  boundary  conditions  must  be 
satisfied.  The  boundary  conditions  used  in  the  present  study  are  listed  in  table  2.  The 
numeral  1  in  this  table  indicates  that  the  nodal  parameter  is  retained  and  0  indicates  that  the 
nodal  parameter  is  set  to  zero. 

For  inclined  (or  curved)  boundaries,  it  is  convenient  to  use  a  modified  set  of  nodal 
parameters  including  normal  and  tangential  components  of  displacements  and  stress  resultants 
at  the  boundary  points,  that  is,  u^^  0^^  where 


r'  ^ 

|«ai 

ju„-i 

iV! 

V  y 

roe/3  i 

Na'i3' 

/ 

[^a/3! 

Ma'/3' 

-  C0S(X^,X^^ 


The  element  equations  at  that  boundary  point  are  modified  accordingly.  For  example,  equa¬ 
tions  (17)  are  modified  as  follows: 

=  P\,  +  Mf.j,  (23 

^ ^ 

where  the  relations  between  and  k|j  are  given  by 

(’-4 


=  llrvrv' 

a  3  a3 


^a',l3'+3  ^^3'  ^a,i3+3 
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where  [k]  ,  [k]  ,  [m]  ,  and  (p}  contain  the  stiffness,  geometric  stiffness,  mass,  and  load 

distributions;  Q  and  contain  the  “generalized”  stiffness  distributions;  and 

are  the  vectors  of  nodal  unknowns  composed  of  the  subvectors  i//j  and  at  the  various 

nodes;  and  the  superscript  T  denotes  transposition.  Note  that  in  the  mixed  models 
(eqs.  (31)),  the  stress  resultants  are  assembled  first. 

The  matrices  and  [s^  are  symmetric,  positive  definite,  and  can  be  banded;  the 

matrices  and  are  banded  symmetric;  and  the  matrix  [0  is  sparse. 

For  stress-analysis  problems,  that  is,  X  =  cc  =  0,  the  governing  equations  of  the  displace¬ 
ment  models  (eqs.  (30))  can  be  solved  by  any  of  the  efficient  direct  techniques  published  in 
the  literature.  (See,  e.g.,  refs.  24  to  26.)  On  the  other  hand,  the  governing  equations  of  the 
mixed  models  can  best  be  solved  by  the  hypermatrix  Gaussian  elimination  scheme.  (See 
ref.  27.) 
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For  eigenvalue  problems,  it  is  convenient  to  modify  the  equations  of  the  mixed  models 
(eqs.  (31))  by  first  eliminating  the  stress  resultants  and  then  rewriting  the  resulting  equations 
in  the  following  form: 


+  CO“ 


(32) 


where 


[s]  =  [sj^  [s]”'  [sj  (33) 

The  matrix  (^sj  is  positive  definite. 

EIGENVALUE  EXTRACTION  TECHNIQUES 

In  the  absence  of  the  external  load  vector  fp)  ,  equations  (30)  and  (31)  define  an  alge¬ 
braic  eigenvalue  problem.  Eor  free-vibration  problems  X  =  (p)  =  0,  the  natural  frequencies 
are  obtained  by  applying  the  subspace  iteration  technique  presented  in  reference  28  to  the 
equations  of  the  displacement  model. 

The  technique  is  based  on  the  use  of  simultaneous  inverse  iteration  with  Gram-Schmidt 
orthogonalization.  The  number  of  vectors  used  in  the  iteration  process  is  more  than  the  eigen¬ 
vectors  required,  but  much  less  than  the  dimensions  of  the  matrices  considered. 

For  the  mixed  models,  the  natural  frequencies  are  obtained  by  applying  the  Sturm 
sequence  technique  with  iterations  to  the  modified  equations  (eqs.  (32)).  In  this  technique 
the  desired  roots  are  first  isolated  by  Sturm  sequence  procedure,  then  the  inverse  iteration  tech¬ 
nique  is  applied  for  the  determination  of  individual  roots  along  with  their  eigenvectors.  (See 
ref.  29.) 

For  bifurcation-buckling  problems,  where  only  the  minimum  buckling  load  parameter  is 
required,  it  is  more  efficient  to  use  the  inverse-power  method  presented  in  reference  30  for 
both  the  displacement  and  mixed  models. 

EVALUATION  OF  STRESS  RESULTANTS 


In  the  mixed  models,  once  the  problem  is  solved,  all  the  stress  resultants  are  readily 
available.  On  the  other  hand,  in  the  displacement  models  the  stress  resultants  are  obtained 
from  the  nodal  displacement  parameters  by  using  the  following  relations: 


9^wV3+p 


(34) 
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(35) 


Qa  =  (SjSWVs  +  wV3+/jl) 

The  stress  resultants  obtained  from  equations  (34)  and  (35)  generally  violate  both  the 
interior  differential  equilibrium  and  the  stress-resultant  boundary  conditions  and  generate  discon¬ 
tinuities  at  the  element  nodes.  Therefore,  in  the  present  study  the  customary  procedure  of 
averaging  contributions  of  contiguous  elements  at  common  nodes  is  followed.  Such  averaging 
is  not  needed  for  the  SQH  element. 

Other  techniques  have  been  suggested  to  improve  the  accuracy  of  the  stress  calculations. 
These  include  the  integral  stress  technique  (ref.  31),  which  is  based  on  least-squares  minimiza¬ 
tion  of  the  stress  error  function  within  each  element,  and  the  conjugate  stress  method  (ref.  32), 
which  uses  biorthogonal  expansion  to  the  displacement  approximation.  Both  these  approaches 
involve  additional  computational  efforts  and  are  not  used  in  the  present  study. 

NUMERICAL  STUDIES 


To  assess  the  relative  merits  of  the  different  displacement  and  mixed  finite-element  mod¬ 
els  developed  in  this  study  (table  1),  a  large  number  of  linear  stress-analysis,  free-vibration,  and 
bifurcation-buckling  problems  are  solved  by  these  finite-element  models.  Particular  emphasis  is 
placed  on  the  effects  of  shear  deformation  and  anisotropic  material  behavior  on  the  accuracy 
and  rate  of  convergence  of  the  different  models. 

The  numerical  examples  are  aimed  at  clarifying  a  number  of  questions  concerning  each 
of  the  following  effects  on  the  accuracy  and  rate  of  convergence  of  finite-element  solutions: 

(a)  an  increase  in  the  order  of  approximating  polynomials,  (b)  addition  of  internal  degrees  of 
freedom,  and  (c)  use  of  derivatives  of  generalized  displacements  as  nodal  parameters. 

PLATE  EVALUATION  RESULTS 

Four  sets  of  plate  problems  are  solved  which  contain  some  of  the  characteristics  typical 
of  practical  problems  and  at  the  same  time  are  problems  for  which  an  essentially  exact  solu¬ 
tion  can  be  obtained.  In  one  of  the  problems,  comparison  is  made  with  experimental  results. 
The  problems  examined  are 

(a)  Stress,  free  vibration,  and  bifurcation  buckling  of  laminated  orthotropic  square  plates 

with  simply  supported  edges 

(b)  Stress  analysis  of  orthotropic  square  plates  with  clamped  edges 

(c)  Stress  and  bifurcation-buckling  analysis  of  square  anisotropic  plates  with  simply  sup¬ 

ported  edges 

(d)  Stress  analysis  of  cantilevered  skew  plates 
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All  the  models  in  table  1  are  applied  to  problems  (a)  and  (b).  The  higher  order  dis¬ 
placement  and  mixed  elements  are  applied  to  problem  (c).  The  higher  order  quadrilateral  dis¬ 
placement  models  SQH  and  SQ12  are  applied  to  problem  (d).  Tlie  results  of  these  studies  are 
discussed  subsequently. 


Square  Plates 

The  first  set  of  problems  considered  is  that  of  the  stress,  free  vibration,  and  bifurcation 
buckling  of  orthotropic  and  anisotropic  square  plates.  Most  of  the  results  presented  in  this 
section  are  for  the  symmetrically  laminated  nine-layered  graphite-epoxy  plates  shown  in  figure  3. 
For  these  plates  two  fiber  orientations  are  analyzed: 

(a)  Orthotropic  plates  with  fiber  orientation  (0/90/0/90/0/90/0/90/0) 

(b)  Anisotropic  plates  with  fiber  orientation  {0 j-O jO j-O jO j-O jO j-O jO),  where  0  <  d  ^  45^ 

For  orthotropic  plates  the  total  thickness  of  the  0^  and  90^  layers  is  the  same,  and  for 
anisotropic  plates  the  total  thickness  of  the  d  and  ~d  layers  is  the  same.  Boundary  conditions 
for  both  simply  supported  and  clamped  plates  are  considered. 

Simply  Supported  Orthotropic  Plates 

The  orthotropic  plate  problems  are  selected  because  an  exact  (analytic)  solution  can  be 
obtained,  and  therefore,  a  reliable  assessment  of  the  accuracy  of  the  different  finite-element 
models  can  be  made.  The  various  solutions  obtained  are  listed  first  and  are  discussed  subse¬ 
quently.  Since  doubly  symmetric  deformations  of  the  plate  are  considered,  only  one-quarter 
of  the  plate  was  analyzed,  and  the  symmetric  boundary  conditions  along  the  center  line  are 
listed  in  table  2. 

For  stress-analysis  problems,  the  plates  were  subjected  to  uniform  loading  p^^.  In  addi¬ 
tion  to  studying  the  accuracy  of  the  maximum  displacements  and  stress  resultants  obtained  by 
the  various  displacement  and  mixed  models,  an  error  index  Ef  (a  function  of  f)  has  been 
introduced  to  provide  a  quantitative  measure  of  the  relative  accuracy  of  the  stress  resultants 
and  displacements  obtained  by  the  different  models.  The  error  index  is  given  by 


Ef  = 


f.  _  f. 
f 

max 


(36) 


where 

f  any  of  the  stress  resultants  or  generalized  displacements 
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f-  f- 

ip  ij 


f 

max 


n 


exact  and  approximate  values,  respectively,  of  the  function  at  the  ith  node 

maximum  absolute  value  of  the  exact  function  in  the  domain  of  interest  (one-quarter 
of  the  plate) 

total  number  of  nodes  in  one-quarter  of  the  plate 


The  error  index  (eq,  (36))  is  essentially  a  weighted  root-mean-square  error.  The  smaller  the 
error  index  E^,  the  more  accurate  the  approximate  solution  (obtained  by  the  finite-element 
model)  is. 

To  study  the  effect  of  shear  deformation  on  the  performance  of  the  different  finite- 
element  models,  three  values  of  the  thickness  ratio  h/a  of  the  plate  were  considered: 
h/a  =  0.1,  0.01,  and  0.001.  As  a  quantitative  measure  of  the  shear  deformation,  the  ratio  of 
the  strain  energy  due  to  transverse  shears  to  the  total  strain  energy  was  computed  for  the  three 
plates.  The  results  are  shown  in  table  3.  As  can  be  seen  from  this  table,  the  shear  deforma¬ 
tion  is  quite  important  for  the  first  plate  and  is  negligible  for  the  latter.  Table  4  gives  the 
values  of  the  error  index  Ef  for  each  of  the  stress  resultants  and  generalized  displacements 
obtained  by  some  of  the  stiffness  and  mixed  finite-element  models  for  two  plate  thicknesses 
(h/a  =  O.l  and  0.01)  and  three  different  grids.  An  indication  of  the  accuracy  and  rate  of  con¬ 
vergence  of  the  solutions  obtained  by  the  different  models  is  given  in  figures  4  and  5,  and  the 
effect  of  h/a  on  the  accuracy  of  the  different  models  is  shown  in  figure  6. 

The  doubly  symmetric  free-vibration  modes  of  the  plate  are  analyzed  by  the  various  ele¬ 
ment  models.  An  indication  of  the  accuracy  and  rate  of  convergence  of  the  fundamental  fre¬ 
quency  obtained  by  different  displacement  and  mixed  models  is  given  in  table  5  and  figure  7 
for  plates  with  thickness  ratios  h/a  of  0.1  and  0.01.  Figure  8  shows  the  effect  of  addition 
of  internal  degrees  of  freedom  on  the  accuracy  and  rate  of  convergence  of  the  four-  and 
eight-node  stiffness  quadrilateral  elements.  Table  6  shows  the  rate  of  convergence  of  the  three 
vibration  frequencies  cuj  3,  j,  and  ^3  3  obtained  by  different  stiffness  models. 

To  study  the  effect  of  the  bending-extensional  coupling  on  the  accuracy  of  the  higher 
order  models,  the  SQ12  and  SQH  elements  were  applied  to  the  free-vibration  problem  of 
two-layered  orthotropic  plates.  Results  obtained  by  these  two  elements  for  the  two  plates 
with  h/a  =  0.1  and  0.01  are  shown  in  table  7  along  with  the  exact  solutions. 

As  a  quantitative  measure  of  the  shear  deformation,  the  exact  frequencies  obtained  by 
the  shear-deformation  and  classical  theories  are  compared  in  tables  5,  6,  and  7. 

Since  the  accuracy  of  the  different  elements  for  buckling  problems  is  expected  to  be 
similar  to  that  for  vibration  problems,  only  the  SQ12  and  SQH  elements  were  applied  to  the 
bifurcation  buckling  of  a  plate  subjected  to  uniaxial  edge  compression  The  results 
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obtained  using  a  2  X  2  grid  in  the  plate  quarter  are  given  in  table  8  along  with  the  exact 
solutions  for  the  three  thickness  ratios  h/a  =  0.1,  0.01,  and  0.001. 

An  examination  of  the  results  obtained  for  simply  supported  orthotropic  plates  reveals 

(1)  Although  the  convergence  of  the  solutions  obtained  by  all  the  displacement  models 
is  monotonic  in  character,  the  convergence  of  the  lower  order  models  is  much  slower  than 
that  of  the  higher  order  models.  This  is  particularly  true  for  stress  resultants  and  for  thinner 
plates.  (See  figs.  4  and  7.) 

(2)  For  the  same  total  number  of  degrees  of  freedom,  the  higher  order  displacement 
models  (e.g.,  SQ12  and  SQH)  lead  to  considerably  more  accurate  results  than  the  lower  order 
models.  This  is  particularly  true  for  stress  resultants  and  for  thinner  plates.  (See  fig.  5.) 

The  same  phenomenon  is  observed  for  vibration  frequencies.  As  an  example  of  this,  for  plates 
with  h/a  =  0.1,  the  fundamental  frequency  obtained  by  the  SQ12  and  SQH  elements  and 
2X2  grid  (corresponding  to  99  and  108  degrees  of  freedom)  agrees  with  the  exact  frequency 
to  four  significant  digits.  (See  table  5.)  In  contrast,  the  error  in  the  fundamental  frequency 
obtained  by  the  SQ4  element  and  5X5  grid  (108  degrees  of  freedom)  is  approximately 
2  percent.  For  higher  frequencies  and  thinner  plates,  the  accuracy  of  the  SQ4  element  dete¬ 
riorated  much  more  rapidly  than  that  of  the  higher  order  models.  (See  tables  5  and  6.) 

(3)  The  accuracy  of  the  solutions  obtained  by  the  lower  order  displacement  models 
(SQ4  element)  is  very  sensitive  to  variations  in  the  thickness  ratio  of  the  plate.  For  thinner 
plates,  the  accuracy  of  this  element  was  found  to  be  very  poor.  (See  tables  4,  5,  and  6.) 

This  is  because  the  assumed  displacement  functions  require  that  the  element  edges  remain 
straight,  and  the  predominant  bending  deformation  in  thin  plates  is  therefore  poorly  represented. 
This  fact  has  been  recognized  by  previous  investigaiors  and  improvements  have  been  suggested. 
(See,  e.g.,  refs.  12,  14,  15,  and  33.)  However,  no  procedure  exists  to  improve  the  accuracy 

of  the  element  for  all  ranges  of  thickness  ratio  of  the  plate. 

(4)  The  SQ8-4  element,  with  different-order  polynomial  approximations  for  displacements 
and  rotations,  although  considerably  more  accurate  than  the  SQ4  element,  is  found  to  be  less 
accurate  than  the  SQ8  element.  (See  fig.  4.)  For  thin  plates  (h/a  =  0.001),  the  performance 
of  the  SQ8-4  element  was  found  to  be  unsatisfactory.  (See  fig.  6.) 

(5)  Of  all  the  finite-element  models  considered,  the  most  accurate  results  for  a  given  total 

number  of  degrees  of  freedom  were  obtained  with  the  SQH  element.  (See  fig.  5  and  tables  5 

and  6.)  The  SQH  element  has  the  added  advantage  that  the  stress  resultants  are  continuous 

along  the  interelement  boundaries  and  no  averaging  is  needed  in  their  evaluation.  However,  in 
the  presence  of  concentrated  loads  or  discontinuities  in  the  geometric  or  material  characteristics, 
some  of  the  nodal  parameters  are  discontinuous  and  a  special  treatment  is  needed.  (See,  e.g., 
ref.  34.) 
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(6)  Bending-extensional  coupling  does  not  appear  to  have  any  adverse  effect  on  the  accu¬ 
racy  of  the  higher  order  displacement  models.  (See  table  7.) 

(7)  The  addition  of  internal  degrees  of  freedom  (bubble  modes)  to  the  displacement  mod¬ 
els  results,  in  general,  in  improving  the  performance  of  the  element.  (See  tables  4,  5,  and  6 
and  fig.  8.)  In  stress-analysis  problems  where  the  internal  degrees  of  freedom  can  be  eliminated 
by  static  condensation  techniques,  this  is  an  effective  way  of  improving  the  accuracy  of  the  ele¬ 
ment,  without  affecting  the  accuracy  of  the  solution.  For  free-vibration  problems,  the  addition 
of  internal  degrees  of  freedom  is  less  effective  than  the  addition  of  nodes  to  the  element.  An 
exception  to  this  is  the  case  of  the  SQ8  element  when  applied  to  the  analysis  of  higher  vibra¬ 
tion  modes  of  plates.  In  this  case  addition  of  higher  order  polynomial  terms  associated  with 
internal  degrees  of  freedom  has  a  more  pronounced  effect  on  the  accuracy  than  the  addition 

of  nodes.  (Compare  the  frequencies  obtained  by  SQ9  and  SQ12  elements  for  the  case  m  =  3, 
n  =  3  in  table  6.) 

(8)  Whereas  for  the  SQ4  element  addition  of  a  single  internal  degree  of  freedom  results 
in  considerable  improvement  in  accuracy,  for  the  SQ8  element  three  internal  degrees  of  freedom 
have  to  be  added  before  a  pronounced  effect  on  accuracy  can  be  observed.  (See  fig.  8.)  An 
exception  to  tliis  is  the  case  of  higher  vibration  modes,  where  the  addition  of  a  single  internal 
degree  of  freedom  improves  the  accuracy  of  the  SQ8  element  substantially.  (See  table  6.) 

(9)  The  solutions  obtained  by  the  mixed  models  are  more  accurate  and  less  sensitive  to 
variations  in  the  thickness  ratio  of  the  plate  than  those  obtained  by  the  displacement  models 
based  on  the  same  shape  functions.  (See  tables  4  and  5  and  figs.  4,  5,  and  6.)  However,  the 
convergence  of  the  solutions  obtained  by  the  lower  order  mixed  models  (MT3  and  MQ4)  is 
slow  and  oscillatory  in  character.  Also,  for  a  given  number  of  degrees  of  freedom,  the  accu¬ 
racy  of  the  solutions  obtained  by  mixed  models  is  lower  than  that  obtained  by  higher  order 
displacement  models  (SQH,  STIO,  and  SQ12).  (See  fig.  5.) 

Two  other  conclusions  were  found  but  the  solutions  on  which  they  are  based  are  not 
reported  herein.  These  are 

(10)  The  accuracy  of  the  solutions  obtained  by  the  triangular  elements  was  found  to  be 
sensitive  to  the  choice  of  their  orientation.  The  best  accuracy  was  obtained  when  the  displace¬ 
ment  models  (ST6  and  STIO)  had  opposite  orientation  to  that  of  the  mixed  models  (MT3 

and  MT6).  (See  fig.  4.)  The  results  shown  in  tables  4,  5,  and  6  and  in  figures  4,  5,  6, 
and  7  were  obtained  for  the  aforementioned  choice. 

(11)  The  effect  of  satisfying  the  force  boundary  conditions  for  the  SQH  element  (in  addi¬ 
tion  to  the  kinematic  conditions),  was  found  to  be  insignificant.  Differences  occurred  only  in 
the  fourth  significant  digit. 
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Before  closing  this  section,  a  comparison  of  the  elements  developed  in  the  present  study 
with  those  previously  reported  in  the  literature  is  in  order.  Since  most  of  the  latter  elements 
do  not  include  shear  deformation,  the  problem  of  an  isotropic  scpiare  plate  with  h/a  =  0.01, 
for  which  the  shear  deformation  is  negligible,  was  selected.  The  plate  had  simply  supported 
edges  and  was  subjected  to  uniform  loading  p^.  The  convergence  of  solutions  obtained  by 
several  classical  plate  elements  was  reported  in  reference  11.  Figure  9(a),  which  is  reproduced 
from  reference  11,  is  contrasted  with  figures  9(b),  (c),  and  (d),  which  show  the  convergence 
of  the  center  displacement  w,  center  bending  moment  Mjj,  and  strain  energy  U  obtained 
by  a  number  of  displacement  and  mixed  shear-flexible  elements.  Except  for  very  coarse  grids 
(2  X  2  or  less  in  the  plate  quarter),  the  higher  order  elements  developed  in  the  present  study 
are  competitive  with  the  refined  elements  previously  reported  in  the  literature.  The  problem 
of  the  thin  isotropic  plate  represents  a  rather  severe  test  for  the  accuracy  of  the  shear-flexible 
elements,  since  the  accuracy  of  such  elements  reduces  with  the  diminishing  of  shear  deformation. 

Clamped  Plates 

To  study  the  effect  of  clamped  edges  as  boundary  conditions  on  the  accuracy  of  the 
different  stiffness  models,  the  edges  of  the  orthotropic  plates  considered  in  the  previous  sub¬ 
section  were  assumed  to  be  totally  clamped  and  the  plates  were  analyzed  by  the  different 
stiffness  and  mixed  models.  The  plates  were  subjected  to  uniform  loading  of  intensity  p^. 

The  standard  of  comparison  was  taken  to  be  the  solution  obtained  by  the  SQH  element  and 
a  6  X  6  grid  in  the  plate  quarter  for  h/a  =  0.1,  and  an  8  X  8  grid  for  h/a  =  0.01 

and  0.001.  An  indication  of  the  accuracy  and  rate  of  convergence  of  displacements  and  stress 

resultants  obtained  by  the  different  models  is  given  in  figure  10  for  three  plate  thicknesses, 
namely,  h/a  =  0.1,  0.01,  and  0.001.  Also,  figure  11  shows  the  distribution  of  the  transverse 

displacement  w  and  the  bending  moment  for  the  thinner  plates  (with  h/a  =  0.01 

and  0.001)  obtained  by  the  higher  order  displacement  models  SQ12  and  SQH  and  the  mixed 
model  MQ8  with  a  2  X  2  grid  in  the  plate  quarter.  As  can  be  seen  from  figure  10,  the 
solutions  obtained  by  the  different  displacement  and  mixed  models  were,  in  general,  less  accu¬ 
rate  than  those  for  simply  supported  edges  (fig.  6).  This  is  particularly  true  for  thinner  plates. 
An  exception  to  this  is  the  SQH  element,  which  exhibited  very  high  accuracy  and  fast  conver¬ 
gence  for  all  thickness  ratios.  Also,  the  remarks  made  in  the  previous  subsection  regarding  the 
effect  of  h/a  on  the  accuracy  and  convergence  of  the  solutions  obtained  by  different  models 
were  found  to  apply  in  this  case,  as  well. 

Anisotropic  Plates 

To  study  the  effect  of  anisotropy  on  the  performance  of  the  higher  order  displacement 
models,  the  fiber  orientations  of  the  graphite-epoxy  plate  shown  in  figure  3  were  chosen  to 
be  {0 j-Q jO  1-6 jO j-6 jd l-d jO)  with  0  <  d  ^  45^.  The  plate  had  simply  supported  edges  and 
was  subjected  to  uniform  loading  of  intensity  p^. 
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Before  the  numerical  studies  were  conducted,  the  effects  of  variations  of  9  on  the 
response  of  the  plate  were  studied.  Also,  an  attempt  was  made  to  introduce  a  quantitative 
measure  of  the  degree  of  anisotropy  of  the  plate.  Since  the  elastic  coefficients 
(with  a  P)  and  ^cxPjp  (with  either  a  =  /3  and  7  p  or 

a  ^  P  and  7  =  p)  vanish  for  orthotropic  (and  isotropic)  plates  and  are  nonzero  only  for 
anisotropic  plates,  it  seems  reasonable  to  take  their  contribution  to  the  total  strain  energy  of 
the  plate  as  a  quantitative  estimate  of  its  degree  of  anisotropy.  Henceforth,  the  contributions 
of  the  anisotropic  coefficients  to  the  total  strain  energy  will  be  referred  to  as  Ug. 

Figure  1 2  shows  the  effect  of  variations  in  9  on  the  values  of  the  displacement  w 

and  the  bending-moment  resultant  Mi  1  at  the  center  of  the  plate  as  well  as  on  the  strain 

^  ^  o 

energies  U,  U^,  and  An  examination  of  figure  12(c)  reveals  that  the  case  6  =  45 

leads  to  the  highest  degree  of  anisotropy  and  the  maximum  value  of  the  shear  deformation. 

Therefore,  the  anisotropic  plate  with  6  =  45^  was  adopted  for  the  convergence  studies. 

An  indication  of  the  accuracy  and  convergence  of  the  higher  order  displacement  mod¬ 
els  STIO,  SQ12,  and  SQH  and  the  mixed  model  MQ8  is  given  in  figure  13  for  the  plate  thick¬ 
nesses  h/a  =  0.1,  0.01,  and  0.001.  The  standard  of  comparison  (converged  solution)  was  taken 
to  be  the  solution  obtained  by  the  SQH  element  and  an  8  X  8  grid  in  the  whole  plate.  Fig¬ 
ure  14  shows  the  distribution  of  the  transverse  displacement  w  and  the  stress  resultant  Mjj 
for  the  thinner  plates  (h/a  =  0.01  and  0.001)  obtained  by  the  SQ12  and  SQH  elements  with 
a  4  X  4  grid,  along  with  the  converged  solutions.  As  in  the  cases  of  simply  supported  and 
clamped  orthotropic  plates,  the  fastest  convergence  was  obtained  by  using  the  SQH  elements. 

The  only  adverse  effect  of  the  anisotropy  on  the  performance  of  the  elements  is  in  the  non¬ 
monotonic  character  of  the  convergence  of  stress  resultants.  (See  fig.  13(b).) 

As  a  further  check  on  the  accuracy  of  the  SQH  elements  in  the  case  of  anisotropic 
plates,  the  bifurcation-buckling  problem  of  the  eight-layered  anisotropic  plate  shown  in  fig¬ 
ure  15  was  analyzed.  The  plate  is  subjected  to  combined  compressive  and  shear  edge  loading. 
The  same  plate  was  analyzed  in  reference  35  using  Galerkin’s  method.  The  results  obtained 
using  three  grid  sizes  of  SQH  elements  (in  the  whole  plate)  are  given  in  table  9  along  with 
those  of  reference  35.  Also,  the  buckling  mode  shapes  are  shown  in  figure  15. 

Skew  Plates 

The  next  problem  considered  is  that  of  the  stress  analysis  of  an  isotropic  skew  plate 
subjected  to  uniform  transverse  loading  (fig.  16).  The  problem  was  selected  because  it  includes 
a  more  complex  set  of  boundary  conditions  and  stress  patterns  than  the  ones  previously 
considered. 

For  this  plate  and  these  boundary  conditions,  an  unbounded  bending  moment  and  a  stress 
singularity  occur  at  point  B.  (See  ref.  36.)  The  nature  of  the  singularity  remains  unaltered 
even  when  the  shear -deformation  theory  (ref.  37)  is  used. 
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Analytical  and  experimental  studies  of  this  problem  were  reported  in  reference  38.  The 
analytic  solution  was  obtained  by  applying  the  mixed  Hellinger-Reissner  formulation  in  conjunc¬ 
tion  with  direct  variational  methods  to  the  classical  plate  theory  (with  shear  deformation 
neglected). 

The  plate  was  analyzed  with  both  the  SQ12  and  SQH  elements.  An  indication  of  the 
accuracy  and  convergence  of  solutions  obtained  by  both  elements  is  given  in  figures  16(a) 
and  (b).  Shown  in  figures  16(c)  and  (d)  are  the  experimental  and  analytical  solutions  of 
reference  38  compared  with  the  present  solutions. 

An  examination  of  figures  16(c)  and  (d)  reveals  that  the  solutions  obtained  by  both  the 
SQH  and  SQ12  elements,  in  addition  to  having  fast  monotonic  convergence,  exhibit  clearly  the 
sharp  gradient  (singularity)  of  the  bending-moment  resultant  M99  at  point  B.  Of  the  two 
finite-element  solutions,  the  SQH  solution  has  a  faster  convergence  and  appears  to  be  more 
accurate.  Moreover,  for  a  4  X  4  or  finer  grid,  the  total  number  of  degrees  of  freedom  in 
the  SQH  solution  is  less  than  those  in  the  corresponding  SQ12  solution. 

SHELL  EVALUATION  RESULTS 

Five  sets  of  shell  problems  are  solved  by  the  displacement  models  developed  in  the  pres¬ 
ent  study.  Comparison  is  made  with  exact  and  other  approximate  solutions  whenever  available. 
These  problems  are 

(a)  Stress  and  free-vibration  analysis  of  orthotropic  shallow  spherical  segments 

(b)  Stress  analysis  of  anisotropic  shallow  spherical  segments 

(c)  Stress  analysis  of  an  isotropic  cylindrical  shell  with  a  circular  cutout 

(d)  Free  vibrations  of  an  orthotropic  cylindrical  shell 

(e)  Free  vibrations  of  an  anisotropic  cylindrical  shell 

All  the  displacement  models  listed  in  table  1  are  applied  to  problem  (a).  Only  the 
higher  order  models  are  applied  to  problem  (b).  The  isoparametric  SQ12  element  is  applied 
to  problem  (c),  and  the  SQH  element  is  applied  to  problems  (d)  and  (e).  The  results  of 
these  studies  are  discussed  subsequently. 

Sliallow  Spherical  Shells 

As  a  first  application  to  a  shallow-shell  problem,  consider  the  stress  and  free-vibration 
analyses  of  simply  supported,  nine-layered,  graphite-epoxy  spherical  segments.  The  geometric 
and  material  characteristics  of  the  shell  are  shown  in  figure  17.  As  for  the  laminated  plates 
examined  in  the  previous  subsections,  shallow  shells  with  two  fiber  orientations  have  been 
analyzed: 
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(a)  Orthotropic  shells  with  fiber  orientation  (0/90/0/90/0/90/0/90/0) 

(b)  Anisotropic  shells  with  fiber  orientation  {6 j-d jO j-d jd j-d jd j-0 jO),  with  0  <  0  ^  45^ 


Orthotropic  Shallow  Shells 

For  the  orthotropic  shells  considered,  analytic  solutions  were  obtained  and  used  as  a 
standard  for  comparing  the  different  finite-element  solutions.  Doubly  symmetric  deformations 
of  the  shell  were  considered,  and  therefore,  only  one-quarter  of  the  shell  was  analyzed. 

For  stress-analysis  problems,  the  shells  were  subjected  to  uniform  loading  p^.  The 
different  displacement  models  were  used  to  obtain  solutions  for  three  thickness  ratios  of  the 
shell  (h/a  =  0.1,  0.01,  and  0.001).  As  a  quantitative  measure  for  the  shear  deformation,  the 
ratios  of  the  strain  energy  due  to  transverse  shear  to  the  total  strain  energy  of  the  shell  were 
computed  for  the  three  shells.  Results  are  given  in  table  10,  and  as  for  orthotropic  plates, 
the  shear  deformation  is  quite  important  for  the  thickest  shell  and  is  negligible  for  the  two 
thinner  shells. 

An  indication  of  the  accuracy  and  rate  of  convergence  of  the  solutions  obtained  by  the 
different  models  is  given  in  figure  18  for  the  shell  with  h/a  =  0.1.  The  effect  of  h/a  on 
the  accuracy  of  the  different  finite-element  solutions  is  shown  in  figure  19.  The  distributions 
of  the  transverse  displacement  w  and  the  stress  resultants  N22  and  Mjj  obtained  by 
the  higher  order  elements  SQ12  and  SQH  with  a  2  X  2  grid  in  the  shell  quarter  are  shown 
in  figure  20  along  with  the  exact  solutions  for  the  two  thinner  shells  (h/a  =  0.01  and  0.001). 

The  first  four  doubly  symmetric  vibration  frequencies  obtained  by  the  different  displace¬ 
ment  models  are  listed  in  table  11  along  with  the  exact  frequencies  for  two  thickness  ratios 
(h/a  =  0.1  and  h/a  =  0.01).  The  solutions  obtained  using  the  SQ4  element  were,  in  general, 
far  removed  from  the  exact  solutions  and  are  not  reported  herein. 

The  orientation  of  the  ST6  and  STIO  elements,  for  optimum  accuracy,  was  found  to  be 
the  same  as  that  for  orthotropic  plate  problems.  (See  fig.  4.) 

An  examination  of  figures  18,  19,  and  20  and  table  11  reveals  that  the  remarks  made 
in  connection  with  the  orthotropic-plate  problems  regarding  the  effectiveness  of  the  higher 
order  models  (STIO,  SQ12,  and  SQH  elements)  and  the  effect  of  internal  degrees  of  freedom, 
apply  in  this  case  as  well.  The  apparent  poor  performance  of  the  different  models  for  the 
case  of  very  thin  shells  (with  h/a  =  0.001)  is  due  to  the  boundary -layer  effects  exhibited  by 
the  stress  resultants  (see  fig.  20),  hence  the  difficulties  (and  nonmonotonicity)  in  convergence 
observed  in  figure  19.  The  convergence  of  the  total  energy  obtained  by  the  higher  order 
models  was  fast  and  monotonic,  even  for  the  very  thin  shell.  (See  fig.  19(d).) 
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Anisotropic  Shallow  Shells 

For  anisotropic  shells  the  fiber  orientations  were  chosen  to  be  {Q j-O jO j-O jO j-Q jd j-d jd) 
with  0  <  0  ^  45^.  The  shells  were  subjected  to  uniform  loading  of  intensity  p^.  The 
quantitative  measures  for  the  degree  of  anisotropy  and  amount  of  shear  deformation  introduced 
for  anisotropic  plates  were  used  for  the  anisotropic  shallow  shells  as  well. 


Figure  21  shows  the  effect  of  variations  in  d  on  the  values  of  the  center  displace¬ 
ment  w  and  the  center  stress  resultants  N99  and  Mjj  for  two  thickness  ratios  of  the 
shell  (h/a  =  0.1  and  0.01).  Also  shown  (fig.  21(d))  are  the  strain  energies  U,  U^,  and 

The  maximum  values  of  and  U^/U  occur  at  different  values  of  0.  This  is  to  be 

contrasted  with  the  anisotropic  plates,  for  which  the  maximum  values  occurred  at  0  =  45^. 

The  accuracy  and  convergence  studies  were  conducted  for  shells  with  0  =  45^.  Fig¬ 
ure  22  gives  an  indication  of  the  accuracy  and  convergence  of  the  center  displacement  w  and 
the  strain  energy  U  obtained  by  the  higher  order  displacement  models  (STIO,  SQ12,  and  SQH) 
for  the  three  thickness  ratios  h/a  =  0.1,  0.01,  and  0.001.  The  standards  of  comparison  (con¬ 
verged  solutions)  were  taken  to  be  the  solutions  obtained  by  the  SQH  elements.  An  8  X  8 
grid  was  used  for  shells  with  h/a  =  0.1  and  0.01,  and  a  10  X  10  grid  was  used  for  shells 
with  h/a  =  0.001.  The  distributions  of  the  normal  displacement  w  and  the  stress  resul¬ 
tants  N92  and  M|j  obtained  by  the  SQ12  and  SQH  elements  with  a  4  X  4  grid  for  the 
thinner  shells  (with  h/a  =  0.01  and  0.001)  are  shown  in  figure  23  along  with  the  converged 
solutions.  As  in  all  the  previous  problems,  the  SQH  solutions  had  the  fastest  convergence. 

The  degradation  of  accuracy  due  to  anisotropy  for  very  thin  shells,  though  not  pronounced  for 
higher  order  displacement  models,  can  be  clearly  seen  by  comparing  the  results  in  figures  20 
and  23. 


Rigid  Body  Modes 

For  shallow  shells,  the  rigid  body  modes  are  trigonometric  in  character  and  therefore 
are  only  approximated  by  the  polynomial  shape  functions  used  in  the  present  study.  To  assess 
the  accuracy  of  the  approximation,  the  eigenvalues  of  the  stiffness  matrices  of  the  various  dis¬ 
placement  models  were  computed  for  the  three  anisotropic  shallow  shells  with  h/a  =  0.1,  0.01, 
and  0.001.  The  lowest  six  eigenvalues  correspond  to  rigid  body  modes;  the  higher  modes  are 
straining  modes.  Table  12  summarizes  the  lowest  seven  eigenvalues,  the  maximum  eigenvalues, 
and  the  traces  of  the  stiffness  matrices  for  the  various  models.  In  all  cases  the  ratio 
was  greater  than  10^,  which  indicates  that  the  rigid  body  modes  are  satisfactorily  represented 
in  these  models. 
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Cylindrical  Shells 


Isotropic  Cylinder  With  a  Circular  Cutout 

Consider  the  stress  analysis  of  an  isotropic  cylindrical  shell  with  a  circular  cutout  sub¬ 
jected  to  a  uniform  axial  tensile  stress  at  its  free  ends.  The  geometric  characteristics  of  the 
shell  and  loading  are  shown  in  figure  24,  The  problem  was  selected  to  assess  the  accuracy  of 
the  isoparametric  SQ12  elements  in  situations  where  high  stress  gradients  and  curved  boundaries 
occur.  The  shell  and  loading  are  doubly  symmetric,  and  therefore,  only  one-quarter  of  the 
shell  was  analyzed. 

An  approximate  analytic  solution  for  the  problem,  assuming  the  cylinder  to  be  of  infi¬ 
nite  length,  was  given  in  reference  39,  where  it  was  shown  that  for  this  shell,  the  shallow- 
shell  approximation  is  valid.  Therefore,  the  use  of  the  SQ12  elements,  with  local  element 
coordinates  coinciding  with  global  shell  coordinates,  is  justified.  A  difference-based  variational 
solution  was  given  in  reference  40,  Finite-element  solutions  using  higher  order  triangular  ele¬ 
ments  were  reported  in  reference  41.  All  the  aforementioned  solutions  were  based  on  the 
classical  shell  theory  (with  shear  deformation  neglected).  Solution  to  a  similar  cylinder  problem 
using  a  refined  grid  of  shear-flexible  quadrilateral  elements  was  reported  in  reference  42. 

Four  graded  networks  with  4  X  4,  5  X  4,  5  X  6,  and  8  X  6  SQ12  elements  were  used 
to  analyze  the  shell.  (See  fig.  25.)  In  an  attempt  to  make  a  rational  choice  for  the  variation 
of  the  grid  size  in  both  the  x^-  and  X2-directions,  a  variable  grid  parameter  f  was  introduced 
(ref.  43  and  fig.  26): 

fr  =  f  (^r+1  - 

where  fj.  is  the  relative  size  of  the  rth  element,  77  refers  to  each  of  the  xj-  and  X2“ 
coordinates,  r?j.  and  77^.+^  are  the  coordinates  of  the  ends  of  the  element,  and  n  is  the 
number  of  elements  in  the  77-direction.  A  second-degree  polynomial  variation  of  fj.  was 
chosen,  that  is, 

fj.  =  a  +  br  +  cr^  (38) 

where  r  is  the  element  number  1  ^  r  ^  n.  The  coefficients  a,  b,  and  c  of  the  poly¬ 
nomial  are  determined  by  specifying  the  relative  sizes  of  the  first  and  last  elements 
and  and  using  the  following  tluee  equations: 

^Err^l-O  (39) 

r=l 
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-  a  +  b  +  c 


(40) 


=  a  +  bn  +  cn“  (41) 

The  characteristics  of  the  grids  used  in  the  present  study  are  shown  in  figure  26. 

The  maximum  stress  concentrations  ^wj^o  strain  energies  obtained  by  the  four 

grids  are  given  in  table  13  along  with  results  of  previous  investigators.  Membrane  stress  dis¬ 
tributions  obtained  by  the  4X4  and  8X6  grids  are  shown  in  figure  27.  The  high  accuracy 
and  rapid  convergence  of  the  solutions  obtained  by  the  isoparametric  SQ12  elements  are  clearly 
demonstrated  by  tliis  example. 


Orthotropic  Cylinders 

The  natural  frequencies  and  mode  shapes  of  orthotropic,  two-layered,  simply  supported 
circular  cylinders  without  axial  restraint  are  studied.  The  problems  are  selected  to  assess  the 
accuracy  of  the  SQH  elements  when  applied  to  laminated  closed  cylinders  with  high  bending- 
extensional  coupling.  The  geometric  characteristics  of  the  shells  studied  are  shown  in  figure  28. 
Shells  with  fiber  orientation  (90/0)  are  analyzed. 

For  these  cylinders  an  analytic  solution  is  obtained  and  is  used  as  a  basis  for  comparison 
of  the  finite-element  solutions.  It  is  found  that  for  this  shell,  the  shallow-shell  (Donnell’s) 
theory  approximation  is  valid.  The  doubly  symmetric  vibration  modes  of  the  cylinders  are 
analyzed  and  the  symmetric  boundary  conditions  along  three  of  the  edges  are  applied.  This 
eliminates  the  axial  rigid  body  mode  of  the  cylinder  and  allows  obtaining  the  vibration  modes 
having  odd  values  of  m  (axial  direction)  and  even  values  of  n  (circumferential  direction). 
Initially  a  uniform  grid  with  2X2  SQH  elements  was  used  to  model  one  octant  of  the  cyl¬ 
inder  (grid  1,  fig.  29);  however,  this  resulted  in  poor  accuracy  for  the  frequencies  and  mode 
shapes  with  n  4.  Subsequently,  the  2X2  grid  was  modified  to  cover  only  one-eighth  of 

the  circumference  (grid  2,  fig.  29).  This  resulted  in  considerable  improvement  in  the  accuracy 
of  the  frequencies  for  n  =  4.  The  frequencies  obtained  by  the  two  grids  are  given  in  table  14 
along  with  the  analytic  solutions  obtained  by  both  the  shear-deformation  and  classical  shallow- 
shell  theories.  This  table  shows  the  decrease  in  accuracy  as  the  element  size-to-wavelength  ratio 
increases  in  the  circumferential  direction,  as  indicated  by  the  increase  of  n.  Numerically,  the 
error  increases  from  less  than  0.5  percent  for  m  =  1,  n  =  2  to  approximately  25  percent 
for  m  =  1,  n  =  4.  The  increased  stiffness  of  the  finite-element  model  due  to  the  larger  ele¬ 
ment  size-to-wavelength  ratio  has  caused  a  greater  increase  in  tlie  error  of  the  finite-element 
analysis  between  the  two  modes.  The  present  example  shows  that  the  SQH  elements  lead  to 
very  accurate  frequencies  provided  the  element  size  is  less  than  half  the  wavelength  of  the 
vibration  mode. 
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Anisotropic  Cylinders 

As  a  final  example,  consider  the  free-vibration  analysis  of  anisotropic  two-layered  circular 
cylinders.  The  shells  have  the  same  characteristics  as  those  for  the  orthotropic  cylinders  dis¬ 
cussed  in  the  preceding  subsection,  except  for  the  fiber  orientation,  which  is  chosen  to  be 
(45/-45). 

Solutions  are  obtained  using  three  grids  with  2  X  4,  4  X  8,  and  6  X  12  SQH  elements 
in  the  whole  cylinder.  (See  fig.  30.)  In  order  to  eliminate  the  axial  rigid  body  mode  of  the 
cylinder,  U]  is  set  equal  to  zero  at  the  center  of  each  grid.  The  fundamental  frequency 
and  associated  mode  shapes  are  shown  in  figure  31.  The  rapid  convergence  of  the  solutions 
obtained  by  the  SQH  elements  is  clearly  demonstrated  by  this  example. 

CONCLUDING  REMARKS 

Several  shear-flexible  finite-element  models  are  applied  to  the  linear  static,  stability,  and 
vibration  problems  of  plates  and  shells.  The  study  is  based  on  the  shallow-shell  theory  with 
effects  of  shear  deformation,  anisotropic  material  behavior,  and  bending-extensional  coupling 
included.  Both  stiffness  (displacement)  and  mixed  finite-element  models  are  considered.  All 
the  elements  examined  are  conforming,  satisfactorily  represent  the  rigid  body  modes,  and 
exhibit  uniform  convergence  for  stress-analysis,  free-vibration,  and  buckling  problems.  Primary 
attention  in  this  study  is  given  to  the  effects  of  shear  deformation  and  anisotropic  material 
behavior  on  the  accuracy  and  convergence  of  different  finite-element  models. 

On  the  basis  of  the  present  study,  the  following  conclusions  seem  to  be  justified: 

1.  Higher  order  displacement  models  (with  cubic  or  bicubic  interpolation  polynomials) 
have  the  following  advantages  over  lower  order  models: 

(a)  The  total  number  of  unknowns  required  for  a  prescribed  level  of  accuracy  is  less 
in  the  higher  order  than  in  the  lower  order  models.  This  is  particularly  true  for  stress 
resultants  and  for  thinner  plates  (with  negligible  shear  deformation). 

(b)  The  performance  of  the  higher  order  models  is  considerably  less  sensitive  to 
variations  in  the  thickness  ratio  and  shear  deformation  than  that  of  the  lower  order 
models. 

2.  The  use  of  derivatives  of  displacements  as  nodal  parameters  (SQH  element)  has  the 
obvious  advantage  that  the  stress  resultants  are  defined  directly  at  the  nodes  and  no  averaging 
is  needed.  In  addition,  this  results  in  improving  the  performance  of  the  element.  However, 
in  the  presence  of  concentrated  loads  or  discontinuities  in  the  geometric  or  elastic  characteris¬ 
tics  of  the  shell,  some  of  the  parameters  will  be  discontinuous  and  a  special  treatment  is 
needed. 
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3.  The  addition  of  internal  degrees  of  freedom  (bubble  modes)  to  displacement  models 
results,  in  most  cases,  in  improving  the  performance  of  the  element.  In  stress-analysis  problems 
where  the  internal  degrees  of  freedom  can  be  eliminated  by  static  condensation  techniques, 

this  is  an  effective  way  of  improving  the  accuracy  of  plate  and  shell  elements  without  affect¬ 
ing  the  accuracy  of  the  solution.  For  free-vibration  (and  buckling)  problems,  tlie  addition  of 
internal  degrees  of  freedom  is  less  effective  than  the  addition  of  nodes  to  the  element.  An 
exception  to  this  is  the  case  of  the  eight-node  quadrilateral  element  when  applied  to  the 
analysis  of  higher  vibration  modes.  In  this  case,  addition  of  internal  degrees  of  freedom  has 
a  much  more  pronounced  effect  on  the  accuracy  than  the  addition  of  nodes. 

4.  If  mixed  models  are  contrasted  with  displacement  models,  the  following  can  be  noted: 

(a)  Tlie  development  of  mixed  models  involves  considerably  less  algebra  than  the 
development  of  displacement  models. 

(b)  The  performance  of  mixed  models  is,  in  general,  insensitive  to  variations  in  the 
thickness  ratio  and  shear  deformation. 

(c)  Use  of  lower  order  interpolation  functions  (linear  or  bilinear)  leads  to  a  medi¬ 
ocre  type  of  performance.  Considerable  improvement  in  the  performance  is  achieved  by 
using  quadratic  shape  functions. 

(d)  For  a  given  number  of  degrees  of  freedom,  the  higher  order  displacement  mod¬ 
els  (with  cubic  or  bicubic  interpolation  polynomials)  lead  to  higher  accuracy  than  the 
mixed  models  with  quadratic  shape  functions.  The  effective  use  of  mixed  models  requires 
the  development  of  efficient  equation-handling  techniques  (e.g.,  based  on  hypermatrix  stor¬ 
age  schemes). 

5.  Whereas  material  anisotropy  was  shown  to  have  an  adverse  effect  on  the  performance 
of  different  displacement  and  mixed  elements,  the  bending-extensionai  coupling  does  not  seem 
to  have  any  pronounced  effect  on  the  accuracy  and  convergence  of  these  elements. 

Langley  Research  Center 

National  Aeronautics  and  Space  Administration 
Hampton,  Va.  23665 
November  10,  1975 
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FUNDAMENTAL  EQUATIONS  OF  SHEAR-DEFORMATION  SHALLOW-SHELL  THEORY 
The  fundamental  equations  of  the  shallow-shell  theory  are  given  in  this  appendix. 
STRAIN-DISPLACEMENT  RELATIONSHIPS 


The  relationships  between  strain  and  displacement  are 
^aP  =  jK^'p  +  ^P^a)  +  ^ap  ^ 

Xap  =  +  ^p^a) 


2ece3  = 

where  are  the  extensional  strains  of  the  reference  surface  of  the  shell;  Xap  ^re  the 

curvature  changes  and  twist;  and  26^,3  are  the  transverse  shearing  strain  components. 

CONSTITUTIVE  RELATIONS  OF  THE  SHELL 

The  relations  between  the  stress  resultants  and  strain  components  of  the  shell  are 

^aj3  ~  ^apyp  ^  ^aPjp  Xyp 

~  ^a|37p  ^yp  ^aPyp  Xyp 


Qa  -  ^a3P3  ^^P3 
The  inverse  relations  are  given  by 

^aP  ~  ^aPyp  '^yp  ^  ^aPyp  ^7p 
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^(3  ®aj37p  ^7p  ^aj^yp  ^7p 


2^03  “  '^a3^3  Q|3 

The  C,  F,  and  D  coefficients  are  shell  stiffnesses  and  the  A,  B,  and 
are  shell  compliances  defined  in  appendix  B. 


coefficients 
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ELASTIC  COEFFICIENTS  OF  LAMINATED  SHELLS 


ELASTIC  STIFFNESSES  OF  THE  LAYERS 

The  nonzero  stiffness  coefficients  and  0^*3133  of  the  kth  orthotropic  layer  of 


the  shell  referred  to  the  directions  of  principal  elasticity  are  given  by 


Xk)/T(k) 


<^1111  "  ^ 


(k)  _  (k)  ^(k)/-(k) 

^1122  ■  ^LT  ^T  /  ^ 


(k) 

^=2222 


(k) 

"=1212 


(k) 

‘^1313 


(k)  _  ^(k) 

'^2323  “  ^TT 

where  the  subscripts  L  and  T  denote  the  direction  of  fibers  and  the  transverse  direction, 
j^Lt  is  Poisson’s  ratio  measuring  the  strain  in  the  T-direction  due  to  a  uniaxial  normal  stress 
in  the  L-direction: 


^TL  ^L  “  '^LT  % 


X  —  1  —  Pi  T  l’ 


LT  ‘'TL 


and  the  superscript  k  refers  to  the  kth  layer. 
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The  stiffness  coefficients  ^a^-yp  ^nd  ‘^a3ii3  satisfy  tlie  following  symmetry  relationships: 
^apyp  ~  ^7po:j3  “  ^j^ajp  ~  ‘^ajSpy 

and 

^a3/33  "  ^P3a3  "  ^30:133  =  ‘-033/3 

If  the  coordinates  are  rotated,  the  elastic  coefficients  c^j^yp  and  trans¬ 

form  as  components  of  fourth-  and  second-order  tensors,  respectively.  Tiie  transformation  law 
of  these  coefficients  is  expressed  as  follows: 

^ap'y'p’  ^  ^aPyp  ^a,a'  ^(3, (S'  ^y,y'  ^p,p' 

and 

‘^a'3/3'3  "  ‘^03/33  ^0,0'  %,/3' 

where  Ca'^'y'p'  ^a'3(S'3  ^he  stiffness-coefficients  referred  to  the  new  coordinate 

system  x^,'  and 

ELASTIC  COEFFICIENTS  OF  THE  SHELL 
The  equivalent  elastic  stiffnesses  of  the  shell  are  given  by 


where  NL  is  the  total  number  of  layers  of  the  shell  and  h|^  and  hj^_j  are  the  distances 
from  the  reference  surface  to  the  top  and  bottom  surfaces  of  the  kth  layer,  respectively.  Tire 
elastic  compliances  of  the  shell  ^a^yp,  ^aPyp^  ^apyp’  ^a3(S3  obtained  by  inver¬ 

sion  of  the  matrix  of  the  elastic  stiffnesses.  (See  ref.  18.) 
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The  shell  stiffnesses  and  compliance  coefficients  satisfy  symmetry  and  transformation  rela¬ 
tions  similar  to  those  of  the  stiffness  coefficients  of  individual  layers. 

The  density  parameters  of  the  shell  are  given  by 


NL  ^h 


(^0-  mi,  =  X)  /  Ps  {}’  ^3’ 
k=l  . 


where  is  the  mass  density  of  the  kth  layer  of  the  shell. 
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SHAPE  FUNCTIONS  USED  IN  PRESENT  STUDY 
QUADRILATERAL  ELEMENTS 

The  expressions  of  the  shape  functions  for  the  different  elements  developed  in  this  study 
in  terms  of  the  quadrilateral  coordinates  (ref.  44)  are  given  in  this  appendix. 

Bilinear  Shape  Functions 

The  shape  functions  for  the  bilinear  approximations  (elements  SQ4  and  MQ4,  see 


(1,-1) 


Sketch  (a) 


where 


(with  a  =  1,2)  are  the  quadrilateral  coordinates  of  node  j. 


Quadratic  Shape  Functions 

The  shape  functions  for  the  quadratic  approximations  (elements  SQ8  and  MQ8,  see 
sketch  (b))  are  given  by 


Corner  nodes 


+ 


Midside  nodes 


a  =  1,3, 5, 7) 


Sketch  (b) 
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Cubic  Shape  Functions 

The  shape  functions  for  the  cubic  approximations  (element  SQ12,  see  sketch  (c))  are 
given  by 

Corner  nodes 

=  ^(i  +  +  hi)  +  h)  -  10)  (j  =  1 ’4,7,10) 

Other  nodes 

^  ~  hi)  ^  ^  2,3,8,9) 

-  ^^1  +  ~  h")  ^  ~  5,6,11,12) 

Sketch  (c) 


Hermitian  Shape  Functions 

The  Hermitian  shape  functions  (element  SQH,  sketch  (d))  used  in  the  present  study  were 
products  of  the  following  set  of  first-order  Hermite  polynomials  (sketch  (e)): 


fi(n  =  -  3?  +  2) 

f2(r)  -  -  r  + 1) 

f3(D  =  -  3?  -  2) 

f4(f)  =  -  1) 


^  f4(?) 
Sketch  (e) 
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9  V  9  V 

If  the  order  of  the  nodal  parameters  at  each  node  is  chosen  to  be  v, 

0  V 

and  ~ — where  v  denotes  any  of  the  fundamental  unknowns,  then  the  shape  func- 
0^1  0^2 

tions  are  given  by 


W'’  -  fi'U  fe(^2 
'  fi+i'Siifs:«2 
=  V'fi  fe+i'i: 

=  fi+i(^l  ff+rh' 

where  the  subscripts  i  and  £  are  functions  of  j  as  follows: 


j 

i 

c  i 

i 

1 

1 

1  i 

5 

3 

1 

9 

3 

3 

13 

1 

3 

Shape  Functions  Associated  With  Nodeless  Variables  (Bubble  Modes) 


Elements  SQ5  and  SQ9 

These  elements  have  one  bubble  mode  given  by 


=  (l  -  ^d)(l  -  ^2^) 


/j  ^5  for  SQ5^ 
\j  =  9  for  SQ9y 
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Elements  SQ7  and  SQll 

These  elements  have  tluee  bubble  modes  given  by 


=  ll(l  -  £l')(l  -  l2^' 


0  =  5  for  SQ7;  j  =  9  for  SQll) 


mJ+2  =  .  _  1^2 


TRIANGULAR  ELEMENTS 

The  expressions  of  the  shape  functions  for  the  different  elements  developed  in  this  study 
in  terms  of  the  triangular  (or  area)  coordinates  given  in  the  following 

sections. 


Linear  Shape  Eunctions 

The  shape  functions  for  the  linear  approximations  (element  MT3,  sketch  (f))  in  terms  of 
triangular  coordinates  are  given  by 


(0,0,1) 


W-*  = 


0  =  1  to  3) 


(1,0,0) 


(0,1,0) 


Sketch  (f) 

Quadratic  Shape  Functions 

The  shape  functions  for  the  quadratic  approximations  (elements  ST6  and  MT6,  sketch  (g)) 
in  triangular  coordinates  are 

Corner  nodes 


-  l)  G  =  2i  -  1;  i  =  1  to  3  and  is  not  summed) 
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Midside  nodes 

(j  =  2i;  i  =  1  to  3  and  is  not  summed;  ^  ~ 


Cubic  Shape  Functions 

The  shape  functions  for  the  cubic  approximations  (element  STIO,  sketch  (h))  in  triangular 
coordinates  are  given  by 


Corner  nodes  (nodes  1 ,  4,  7) 


A/J  =  ^(3^i  -  l)(3?i  -  2)^i 

(j  =  3i  ~  2;  i  =  1  to  3  and  is 

not  summed) 

®c 

(D 

(D  (3)  0 

Sketch  (h) 

Boundary  nodes 

n'  =  1  Si  £i+i(3Si  -  0  (j 

i  =  3i  -  1;  nodes  2,  5,  8)  , 

(i  =  1  to  3 

and  is  not 

(  summed; 

i 

U  =  li) 

w'  =  f  £i  £i+i(3£i+i  -  1)  (i 

=  3i;  nodes  3,  6,  9)  i 

J 

Interior  node  (node  10) 


M-i  =  275,^2^3 
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FORMULAS  FOR  COEFFICIENTS  IN  GOVERNING 
EQUATIONS  FOR  INDIVIDUAL  ELEMENTS 

The  expressions  for  the  independent  stiffness  coefficients  in  equations  (17)  are  given  by 

■  X(e)  ““ 

K»3  -  /(e)  Wo  Hl> 


=  L,,  Wo  ’oN' 


K33  =  /(e)  if<«rPo  “or  ‘‘Jo  W3 


Ki3,3  =  /(e,  CWo  ‘‘Jo  *  W3 


C3JJ+3  -  /(e)  CWo  “t"*  ®0«'  W3 


The  independent  nonzero  geometric  stiffness  coefficients  are  given  by 


K33  =  /,e)  nSj  3jN' 


The  independent  nonzero  consistent  mass  coefficients  are  given  by 


Ma/3  ■  "^0  ^a|3 


M(i,(3+3 


/■ 

=  /  mi 


5^/.  W*^  dr2 


W33  ==  /"(e)  *"0 
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Ma+3,/3+3  =  >^2 

where  5^  is  the  Kronecker  delta  on  a  and  jS. 

The  expressions  for  the  “generalized”  stiffness  coefficients  in  equations  (18)  are  given  by 

^a+|3-l,7+p-l  “  ^(e)  ^‘^yP  ^ 

^a+(3-l,7+p+2  ~  J'^{e)  ^‘^^yP  ^ 

^a+l3+2,7+p+2  “  ^(e)  ^^&yP  ^ 

^i+6,|3+6  "  J^(g)  ^0i3P3 
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In  the  above  equations  the  contributions  of  the  line  integrals  have  been  neglected  for 
simplicity;  k  is  a  constant  equal  to  1  when  a  and  1/2  when  a  -  p;  the  range  of 

the  lowercase  Latin  indices  is  1  to  m,  where  m  is  the  number  of  shape  functions;  the 
range  of  the  Greek  indices  is  1,2;  and  a  repeated  index  denotes  summation  over  the  full 
range  of  the  index. 

It  should  be  mentioned  that  for  elements  with  internal  degrees  of  freedom  (SQ5,  SQ7, 
SQ9,  and  SQll),  the  indices  ij  in  the  expressions  for  P^  and  P^  were  assumed  to  have 
a  range  equal  to  the  number  of  nodes  in  the  element  (i.e.,  4  for  SQ5  and  SQ7  elements, 
and  8  for  SQ9  and  SQll  elements).  This  means  that  the  loading  was  distributed  on  the  nodes 
of  these  elements  and  no  loading  was  associated  with  internal  degrees  of  freedom. 
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TABLE  1.-  CHARACTERISTICS  OF  SHEAR-FLEXIBLE  FINITE-ELEMENT  MODELS 


Degrees  of  freedom 


TABLE  2.-  BOUNDARY  CONDITIONS  USED  IN  PRESENT  STUDY 
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TABLE  3.-  EFFECT  OF  THICKNESS  RATIO  h/a  ON  TOTAL 
AND  TRANSVERSE  SHEAR-STRAIN  ENERGIES  OF  PLATES 

Simply  supported,  nine-layered,  square  orthotropic  plate  subjected  to 
uniform  pressure  loading  p^;  U  denotes  total  strain  energy  of 
_  plate  and  Ugj.^  denotes  shear-strain  energy  of  plate 
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TABLE  6.-  CONVERGENCE  OF  NONDIMENSIONAL  FREQUENCIES 
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Numbers  in  parentheses  refer  to  classical-theory  solutions  (with  both  shear  deformation  and  rotary  inertia 


TABLE  6.-  Concluded 


TABLE  7.-  ACCURACY  OF  VIBRATION  FREQUENCIES  OBTAINED 
BY  SQI2  AND  SQH  ELEMENTS 

Simply  supported,  two-layered  orthotropic  plate;  2X2  grid; 

V,n  " 


Values  of 

[>m,n  (a/h)0  X  10"*  for  - 

m,n 

SQ12 

SQH 

Analytic  solution 

_ (a) _ 

h/a  = 

0.1 

1,1 

1.0580 

1.0578 

1.0578 

(1.1244) 

1,3 

3,1 

j  4.8447 

4.8405 

4.8305 

(6.5730) 

3,3 

7.3000 

6.8895 

6.8757 

(9.6664) 

h/a  = 

0.01 

1,1 

1.1463 

1.1303 

1.1300 

(1.1308) 

1,3 

3,1 

j  7.1596 

6.9498 

6.7319 

(6.7644) 

3,3 

15.7435 

10.4152 

10.1118 

(10.1725) 

^Numbers  in  parentheses  refer  to  classical-theory 
solutions  (with  both  shear  deformation  and  rotary 
inertia  neglected). 


TABLE  8.-  ACCURACY  OF  BUCKLING  LOAD  PARAMETER  X 
OBTAINED  BY  SQH  AND  SQ12  ELEMENTS 


Simply  supported,  nine-layered,  square  orthotropic  plate  subjected 
to  uniaxial  edge  compression;  N^j  =  -1;  2  X  2  grid  in  one- 
^  quarter  of  plate 


Values  of  Xa^y 

for  - 

h 

a 

SQH 

SQI2 

Analytic  solution 

(a) 

0.1 

27.012 

27.014 

27.0069  (36.1597) 

.01 

36.051 

36.419 

36.0365  (36.1597) 

.001 

36.177 

i 

69.060 

36.1585  (36.1597) 

^Numbers  in  parentheses  refer  to  classical-theory 


solutions  (with  shear  deformation  neglected). 


TABLE  9.-  CONVERGENCE  OF  BUCKLING  LOAD  PARAMETER  X 
OBTAINED  BY  SQH  ELEMENTS 

c~  ^ 

Simply  supported,  eight-layered,  anisotropic  plate  with  fiber 

orientation  (90/0/-45/45/45/-45/0/90)  subjected  to 

combined  compressive  and  shear  edge  loadings; 

^  =  0.0072;  N^j  =  0;  N^^  =  L  ^22  " 


Grid  size 
(full 
plate) 

Values  of  Xa^/Ejh^  for  - 

SQH  element 

Galerkin’s 
method  (ref.  35)® 

2X2 

19.745 

19.590 

3X3 

19.194 

4X4 

19.046 

^Based  on  classical  theory  (with  shear  deformation 
neglected). 


TABLE  10.-  EFFECT  OF  THICKNESS  RATIO  h/a  ON  TOTAL  AND  TRANSVERSE 

SHEAR- STRAIN  ENERGIES  OF  SHELLS 

"simply  supported,  nine-layered,  orthotropic  shallow  spherical  shells 
(R/a  =  10,  f/a  =  0.0125)  subjected  to  uniform  pressure  loading 
p^;  U  denotes  total  strain  energy  and  denotes  shear-strain 

_  energy  - 
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TABLE  11.-  CONVERGENCE  OF  NONDIMENSIONAL  FREQUENCIES 
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Numbers  in  parentheses  refer  to  classical-theory  solution  (with  shear  deformation  and  rotary  inertia 


TABLE  12.-  EIGENVALUES  B  OF  THE  STIFFNESS  MATRICES  FOR  VARIOUS  DISPLACEMENT  MODELS 


TABLE  13.-  COMPARISON  OF  SOLUTIONS  OBTAINED  BY  ISOPARAMETRIC  SQ12 
ELEMENTS  WITH  THOSE  OF  PREVIOUS  INVESTIGATORS 


[circular  cylindrical  shell  with  a  circular  cutout  loaded  in  tensio^ 


Number  of 

‘^ll 

/"o 

_ 

UE/(ao\3) 

(b) 

_4 

Grid 

degrees  of 
freedom 

Membrane 

Membrane  + 
bending 

X  10 

4X4 

525 

3.643 

4.268 

4.734951 

5X4 

640 

3.691 

4.252 

4.735373 

5X6 

920 

3.712 

4.257 

4.711623 

8X6 

1415 

3.666 

4.223 

4.724903 

Finite  differences 
(ref.  40)^ 

753 

3.603 

4.096 

Finite  elements 
(ref.  41)^ 

367 

3.690 

4.249 

4.729269 

Analytic  solution 
(ref.  39)2 

3.658 

4.180 

^Based  on  the  classical  theory  (with  shear  deformation  neglected). 
'^Strain  energy  in  one-quarter  of  the  shell. 


TABLE  14.-  ACCURACY  OF  VIBRATION  FREQUENCIES  OBTAINED 

BY  SQH  ELEMENTS 

Simply  supported,  two-layered,  orthotropic  circular  cylinder 
with  h/R  =  0.05,  R/L|  =  0.5;  \ri,n  ^ 


m,n 

Values  of 

^m,n  for- 

SQH  element 

Analytic 

solution 

(c) 

Grid  1^ 

Grid  2*^ 

1,2 

0.5512 

— 

0.5487  (0.5494) 

1,4 

.7932 

0.6396 

.6356  (  .6473) 

3,2 

1.7173 

- - 

1.7121  (1.7237) 

3,4 

1.4143 

1.3390 

1.3317  (1.3581) 

^Grid  1 :  2  X  2  in  shell  octant. 

^Grid  2:  2  X  2  (2  elements  in  one-eighth  of 
the  circumference). 

^Numbers  in  parentheses  refer  to  classical-theory 
solution  (with  both  shear  deformation  and  rotary  inertia 
neglected). 
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Number  of  degrees  of  freedom  (quarter  plate) 

Figure  5.-  Convergence  of  w  and  obtained  by  different  stiffness  and  mixed 

models  with  increasing  number  of  degrees  of  freedom.  Simply  supported,  nine¬ 
layered,  orthotropic  square  plate  with  ^  =  0.1  and  0.01. 
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refinement.  Simply  supported,  nine-layered,  orthotropic  square  plate. 


Number  of  elements  per  side  (quarter  plate) 

(a)  ^  =  0.1;  four-node  quadrilaterals, 
a 


Figure  8.-  Effect  of  internal  degrees  of  freedom  on  accuracy  and  convergence  of  vibration 
frequencies.  Simply  supported,  nine-layered,  orthotropic  square  plate. 


Number  of  elements  per  side  (quarter  plate) 


(c)  -  =  0.01;  eight-node  quadrilaterals. 
Figure  8.-  Concluded. 


Number  of  elements  per  side,  n  (whole  plate) 


(a)  w  at  center  for  previously  developed  (b)  w  at  center  for  present 

classical-theory  elements  (ref.  11).  shear -flexible  elements. 


Number  of  elements  per  side,  n  (whole  plate) 

(c)  at  center  for  present  elements.  (d)  Strain  energy  for  present 

elements. 

Figure  9.-  Convergence  of  w,  and  U  with  grid  size  for  shear -flexible  elements 

of  present  study  and  some  previously  developed  classical-theory  elements  (ref.  11). 
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;  10.-  Effect  of  h/a  on  convergence  of  bending-moment  resultant  and  transverse  displacement 

obtained  by  different  stiffness  and  mixed  models.  Clamped,  nine-layered,  orthotropic  square  plate. 


Number  of  elements  per  side  (quarter  plate) 

Maximum  bending-moment  resultant  at  edge 

Figure  10.-  Concluded. 


plate  subjected  to  uniform  normal  pressure  loading. 


a  function  of  fiber  orientation. 


h/a  Symbol  x 


3.7323 


2.4079 


Number  of  elements  per  side  (whole  plate) 

(a)  Transverse  displacement  w  at  center. 

Figure  13.-  Effect  of  h/a  on  convergence  of  bending-moment  resultant 

transverse  displacement  w,  and  strain  energy  U  obtained  by  higher  order 
stiffness  and  mixed  models.  Simply  supported,  nine-layered,  anisotropic 
square  plate  (45/-45/45/-45/45/-45/45/-45/45). 


SQH 


Figure  14.-  Distribution  of  transverse  displacement  w  and  bending-moment 
resultant  along  X2  =  Simply  supported,  nine-layered,  anisotropic 

square  plate.  =  0.01  and  0.001. 


Contour  plot  for  <pi/(p^ 


Buckling  mode  shapes  for  laminated  anisotropic  plate  subjected  to  combined  compressive  and 

shear  edge  loadings. 


xperi merit  (ref.  38) 


-  Convergence  of  maximum  transverse  displacement  w  and  bending-moment  resultant  M22  obtained 
by  higher  order  stiffness  models.  Cantilevered  isotropic  skew  plate  under  uniform  loading. 


16.-  Concluded, 


^3 


=  0.25 

Figure  17.-  Characteristics  of  laminated  graphite-epoxy  shallow  shells 

used  in  present  study. 
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Number  of  elements  per  side  (quarter  shell) 


(a)  Transverse  displacement  w  at  center. 

Figure  19.-  Effect  of  h/a  on  convergence  of  transverse  displacement,  stress 
resultants,  and  strain  energy  obtained  by  different  stiffness  models.  Simply 
supported,  nine-layered,  orthotropic  shallow-spherical  segment. 


Number  of  elements  per  side  (quarter  shell) 

(b)  Membrane  stress  resultant  N22  at  center 
Figure  19.-  Continued. 


Number  of  elements  per  side  (quarter  shell) 

(c)  Bendir^ -moment  resultant  at  center 

Figure  19,-  Continued. 
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21.“  Concluded. 


I 


by  different  stiffness  models.  Simply  supported,  nine-layered,  anisotropic  spheric 
segment  with  fiber  orientation  (45/-45/45/-45/45/-45/45/-45/45). 
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Figure  24.-  Cylinder  with  a  circular  cutout  loaded  in  tension 


5x6  Grid  8x6  Grid 

Figure  25.-  Grids  used  in  present  study  for  cylinder  with  a  cutout 


Figure  29.-  Grids  and  modes  for  orthotropic  cylinder 


for  anisotropic  cylinder. 


